Body mass, temperature, and depth shape the maximum intrinsic rate of population increase in sharks and rays

Abstract An important challenge in ecology is to understand variation in species' maximum intrinsic rate of population increase, r max , not least because r max underpins our understanding of the limits of fishing, recovery potential, and ultimately extinction risk. Across many vertebrate species, terrestrial and aquatic, body mass and environmental temperature are important correlates of r max . In sharks and rays, specifically, r max is known to be lower in larger species, but also in deep sea ones. We use an information‐theoretic approach that accounts for phylogenetic relatedness to evaluate the relative importance of body mass, temperature, and depth on r max . We show that both temperature and depth have separate effects on shark and ray r max estimates, such that species living in deeper waters have lower r max . Furthermore, temperature also correlates with changes in the mass scaling coefficient, suggesting that as body size increases, decreases in r max are much steeper for species in warmer waters. These findings suggest that there are (as‐yet understood) depth‐related processes that limit the maximum rate at which populations can grow in deep‐sea sharks and rays. While the deep ocean is associated with colder temperatures, other factors that are independent of temperature, such as food availability and physiological constraints, may influence the low r max observed in deep‐sea sharks and rays. Our study lays the foundation for predicting the intrinsic limit of fishing, recovery potential, and extinction risk species based on easily accessible environmental information such as temperature and depth, particularly for data‐poor species.


| INTRODUC TI ON
Ecologists often have to inform policy decisions with knowledge and experience rather than data. A classic and ongoing challenge is that of making a decision about whether a species is overfished or not. Overexploitation is the leading cause of extinction risk in the ocean and a major cause also on land (Maxwell et al., 2016;Reynolds et al., 2005). A large body of work shows that maximum size, either length or weight is a simple heuristic that can guide our thinking about extinction risk. Setting aside the desirability or catchability of different-sized fishes, broadly speaking, the larger-bodied species have slower life histories and tend to decline faster when fished (Jennings et al., 1998;Reynolds et al., 2005). Body size data are among the most widely available traits, but it is only one dimension of the life histories of indeterminate-growing species (Chichorro et al., 2019). The second dimension of life histories tends to be comprised of time-related "speed of life" traits, such as age at maturity, maximum age, and von Bertalanffy growth rates (Juan-Jordá et al., 2013). In the few cases where high-quality data can be found, it is increasingly clear that speed of life traits are the better correlate of population dynamics and extinction risk (Anderson et al., 2011;Juan-Jordá et al., 2015).
The maximum intrinsic rate of population increase (of species at small population sizes) r max is the integration of speed of life traits and "is perhaps the most fundamental parameter in population biology" (Myers & Worm, 2005). Estimates of r max are commonly used limit reference points as they are equivalent to the level of fishing mortality that will drive a species to extinction and defining the maximum rate of population recovery (Dulvy et al., 2004;Hutchings et al., 2012;Myers et al., 1997). Ever larger databases of life-history traits make it easier to calculate r max , providing an opportunity to seek ecological explanations for variation in the maximum intrinsic rate of population increase. Indeed, Southwood (1977) viewed habitat quality as the template shaping population growth rates. In the ocean, the scope for growth changes profoundly with temperature, such that tropical species have faster life histories and dynamics (Juan-Jordá et al., 2013;Munch & Salinas, 2009) but attain smaller sizes (Cheung et al., 2013;Fisher et al., 2010;Jennings et al., 2008).
From a metabolic perspective alone, at high temperature, fishes are squeezed between low oxygen solubility and high oxygen demand due to temperature-forced metabolic rates (Pauly, 2021;Pörtner et al., 2017;Rubalcaba et al., 2020). Furthermore, within a given thermal habitat, aquatic species exhibit a range of oxygen demands as a result of differing lifestyles, activity patterns, and trophic levels, which may have consequences for life-history and population dynamics (Killen et al., 2010(Killen et al., , 2016Wong et al., 2021). Similar patterns have been observed on land, where allometric variation in mammal production rates (analogous to r max ) can also be shaped by differences in lifestyle and trophic level (Sibly & Brown, 2007, note that temperature is not an axis of variation in endotherms).
The sharks, rays, and chimeras (class Chondrichthyes, hereafter referred to as "sharks and rays") are an ideal taxon for studying the aquatic allometry of demographic rates as they encompass a broad range of sizes and occur across a wide range of temperatures and habitats. The lack of a pelagic larval stage in sharks and rays allows for a more straightforward estimation of population productivity with limited life-history data than their bony counterparts (Myers & Mertz, 1998;Pardo et al., 2016). The maximum intrinsic rate of population increase r max among sharks and rays is known to decrease with increasing size (Dulvy, Pardo, et al., 2014;Hutchings et al., 2012) and depths (García et al., 2008;Simpfendorfer & Kyne, 2009). However, the previous work has only considered categorical habitat classifications, and there is some evidence that the relationship between r max and body size breaks down in the deep sea as even small deep-water sharks and rays have very low r max (Forrest & Walters, 2009;Rigby & Simpfendorfer, 2015), suggesting there are stronger constraints on the mass scaling of r max occurring at greater depths than those imposed by body mass and temperature. Thus, depth provides another environmental correlate of interspecific variation in r max .
Here, we examine the separate roles of temperature and depth on the maximum intrinsic rate of population increase and mass scaling relationship among sharks and rays by using estimates of temperature-at-depth derived from species distributions maps (thermal habitat templates) with an information-theoretic approach, while accounting for phylogenetic nonindependence.

| Data
To determine the role of temperature, depth, and the scaling of body mass on maximum intrinsic rate of population increase, r max of marine sharks and rays, we gathered three types of data: life-history parameters (maximum weight and r max ), species-specific environmental traits (median depth and temperature of occurrence), and phylogenetic trees to account for evolutionary relationships in the model-fitting process. Maximum reported body mass (in grams) was obtained from FishBase (Froese & Pauly, 2016) using the rfishbase package (Boettiger et al., 2012). We used maximum weight instead of weight at maturity as it is more readily available in the literature, and at first approximation, they are proportional to each other. When body mass data were unavailable, maximum length data were converted to body mass using species-specific lengthto-weight conversions also sourced from FishBase. Data for r max were obtained from a modified Euler-Lotka model following Pardo et al. (2016). Median depth estimates for each species as reported by Dulvy, Fowler, et al. (2014), and temperature-at-depth was derived from species distribution maps. These maps were obtained from AquaMaps, which is an online resource of global species distribution models for over 25,000 aquatic species (Kaschner et al., 2015).
The core distributions for each species (i.e., where probability of occurrence ≥0.9) were overlaid with the International Pacific Research Center's interpolated dataset of gridded mean annual ocean temperatures across 27 depth levels (0-2000 m below sea level), which is based on measurements from the Argo Project (see http://apdrc.soest.hawaii.edu/proje cts/Argo/data/stati stics/ On_stand ard_level s/Ensem ble_mean/1x1/m00/index.html for data and more information). This temperature interpolation covers most of the world's oceans but has incomplete coverage of some shallow coastal areas, such as the Indo-Pacific Triangle and the southeast coast of South America.
To calculate a temperature that characterizes the thermal distribution of each species, we selected the depth level from the grid that was closest to the species' median depth, and from that grid extracted all the temperature grid points that overlaid the species' core distribution. From this distribution of temperature values, we calculated the median and set as the temperature value for each species. For example, if a species' median depth was 130 m, we used temperatures from the 150 m depth layer to estimate the median annual temperature as 130 m is closer to 150 than 100 m. In species that are known to be mesotherms (e.g., family Lamnidae), we added a correction factor of 3.5°C.
The phylogenetic tree was obtained from Stein et al. (Stein et al., 2018) and we followed their scientific nomenclature.

| Metabolic scaling expectations
The availability of consistent estimates of environmental temperature derived from species distribution maps allows us to draw from metabolic scaling theory. Savage et al. (2004, see equation 8) showed that across multiple taxa, r max is related to body mass and temperature, where r max is the maximum intrinsic rate of population increase (in year −1 ), i 0 is a taxon-specific normalization constant, M is the adult body mass of a species (in grams), β is the scaling exponent, E is the activation energy, k B is the Boltzmann constant (8.617 × 10 −5 eV), and T is temperature (in Kelvin). The equation above can be simplified when transformed to log space: which is equivalent to a simple linear model where the intercept β 0 is the log-transformed normalization constant, the coefficient β 1 is the mass scaling, and β 2 is the negative activation energy -E, which is the coefficient of inverse temperature 1∕k B T. We do not use a taxon-specific normalization constant (i 0 = β 0 ) as we are using a phylogenetic covariance matrix instead of a taxonomic nested structure. Equation 3 is what would be expected from metabolic scaling theory, and is one of the multiple hypotheses we compare in this study, which are laid out in Table 1.

| What is the role of temperature and depth in the mass scaling relationship of r max ?
We compared nine different models, each representing a different hypothesis of how r max might vary with body mass, depth, and temperature, using an information-theoretic approach (Burnham & Anderson, 2002). We scaled and centered the data for temperature and depth. (1) TA B L E 1 Models of correlates of the maximum intrinsic rate of population increase (r max ) tested in our analysis and the hypotheses associated with each Note: All models were implemented using a phylogenetic generalized least squares framework, where the phylogeny is incorporated as a random effect in the form of a covariance matrix.
To assess whether the estimated median temperature for each species was indeed informative, we also included a model where r max is explained only by variation in body mass. We compared all models using the corrected Akaike Information Criterion (AICc). If including a parameter improved a model's AICc value by <2, we considered that parameter to be uninformative (Arnold, 2010).
To account for nonindependence among closely related species, we included a phylogenetic covariance matrix as a random effect in a generalized least squares (GLS) framework using the pgls function in the caper package, (Freckleton et al., 2002;Orme et al., 2013), running on R version 3.3.2 (R Core Team, 2016). The phylogenetic covariance matrix, which is contrasted with the residuals from each model, was adjusted using Pagel's λ, which is a scaling parameter (i.e., a multiplier of the off-diagonal elements of the covariance matrix), that is estimated in the modeling approach and transforms the phylogenetic covariance matrix so that it best fits the trait data (Freckleton et al., 2002;Pagel, 1999). This allows for a measure of phylogenetic signal ranging from no phylogenetic signal at all (i.e., classic OLS; λ = 0) to the traits evolving under the Brownian motion model (λ = 1) (Revell, 2010).
We tested for collinearity by estimating variance-inflation factors (VIF) for all coefficients in the models explored using the car package (Fox & Weisberg, 2011). None of the VIF values estimated were greater than 2 (except when interactions are included, which is expected). However, depth and temperature are positively correlated (Pearson's r = .64); although this correlation is lower than the threshold of |r| > .7, where collinearity severely distorts model estimation (Dormann et al., 2013).

| RE SULTS
We collected temperature, body mass, depth, and r max data for 63 chondrichthyan species including 40 sharks, 20 rays, and three chimeras ( Figure 1).
The model with the greatest support included temperature and depth as well as an interaction term between body mass and temperature (Table 2). This model also explained the greatest amount of variation (adjusted R 2 = .33).
The model that included the same variables as the top-ranked model as well as an interaction term between body size and depth, received one-third the support of the top-ranked model. While there was partial evidence of a negative interaction between body mass and depth, this model had a ΔAICc value greater than 2 without an increase in adjusted R 2 . Therefore, we considered the interaction between body size and depth to be uninformative and the model with this interaction was not considered further (Arnold, 2010).
while there was some evidence of a negative relationship for some species, this relationship wasn't supported for the majority, hence you proceed using the top-ranked model while recognizing it does not cover all cases.
All three other models that had marginal support (5 ≤ Δ AICc ≤ 7) also included temperature: one was similar to the top-ranked model but lacked the effect of depth; one included a main effect of temperature, while the effect of body mass varied with depth; and the other includes all the main effects with no interactions ( Table 2).
The effect size of body mass in all models is around −0.3 (see shaded areas in Figure 2, Table 3), which overlaps with the expectation of −1∕4 predicted by metabolic scaling theory . Intercepts also consistently overlap zero for all models, albeit with high uncertainty (Table 3). The inclusion of only an additive effect of temperature to the relationship between r max and body mass was uninformative (ΔAICc10 → 9.8), but the inclusion of the interaction between body mass and temperature resulted in an improved model fit (ΔAICc = 5.0; Table 2). By contrast, the addition of additive depth alone or an interaction with depth did not improve the relationship between r max and mass (ΔAICc11 → 9.2); however, the addition of depth to the model with an interaction between body mass and temperature improved this model's support twelvefold (weight of 0.63 vs. 0.052, Table 2). This effect of depth is consistently negative across the top-ranked models, indicating that r max decreases with increasing depth.
The interaction between body mass and temperature was always positive (Table 3)  There was a strong phylogenetic signal from the residuals of r max in all ten models examined. Estimated λ ranged between 0.78 and 0.87 with the best model having the lowest λ value (Table 3).
The phylogeny by Stein et al. (2018) did not provide a single" best" tree but a distribution of possible trees with the same topology but differing branch lengths. We ran our analyses by using 20 different trees sampled from their distribution. The results were almost identical regardless of which tree was used, and therefore we only report our findings when using only a single tree (Table 4).

| DISCUSS ION
We generally find that the maximum intrinsic rate of population increase is lower in larger-bodied species. However, we also found two intriguing depth-and temperature-related effects. First, we found that r max decreases at greater depths across marine sharks and rays.
Further, the relationship between mass and r max is temperaturedependent, such that r max decreases with increasing body mass more steeply at higher temperatures. We also found that the scaling of body mass and r max changes across a temperature gradient, which is in line with the observation by Rigby and Simpfendorfer (2015) and suggests the absence of a relationship between size and r max in deep-sea sharks and rays is perhaps due to the low temperatures found in the deep sea rather than depth itself. In other words, it is temperature rather than depth that "flattens" the mass scaling of r max in sharks and rays, and species that live in colder environments have a shallower mass scaling relationship with r max than those that live in warmer waters. Overall, the interaction between the effect of body mass and temperature likely drives the nonlinear scaling of shark and ray r max . Our findings imply that body size might not be a good indicator of r max among cold-water shark and ray species as their body mass coefficient approaches zero, which is an important consideration for the conservation and management of these often data-poor species.
There is increasing evidence that the allometries of biological rates can be nonlinear, which can stem from multiple separate mechanisms acting concurrently (Glazier, 2005(Glazier, , 2015Kolokotrones et al., 2010). Thus, differences in temperature could potentially affect scaling relationships via separate mechanisms (Bruno et al., 2015;Ohlberger et al., 2012). The differing effects of physiological, energetic, and geometric processes on biological rates across allometries could explain the nonlinearity of the relationship observed in this study and are discussed below.
As dissolved oxygen concentrations decrease as water temperature increases (Pauly, 2010), the largest bodied species and F I G U R E 1 Phylogeny, body mass, maximum intrinsic rate of population increase (r max ), temperature, and median depth in sharks, rays, and chimeras. Phylogenetic tree is based on Stein et al. (2018), body mass estimates were sourced from FishBase (Froese & Pauly, 2016), r max estimates from Pardo et al. (2016), median depth values from Dulvy, Fowler, et al. (2014), and mean annual temperature values (at median depth) were estimated based on species distribution maps from AquaMaps (Kaschner et al., 2015) and global temperature grids from the Argo database. Asterisks (*) denote temperature values corrected for mesothermy. Horizontal dotted lines indicate separate taxonomic orders.
individuals may become more oxygen limited, which potentially explains why r max decreases more steeply with increasing body size in warmer water species than colder water ones (Rubalcaba et al., 2020). One mechanism that might plausibly account for this is the gill-oxygen limitation theory (GOLT, Pauly & Cheung, 2018;Pauly, 2021). According to this theory, gill surface area imposes a stronger limitation on oxygen uptake as organisms grow larger, which results in less energy allocated to growth and reproduction in warmer and also deeper waters (i.e., waters with less dissolved oxygen). While physiologists dispute the directionality of this mechanism (e.g., Lefevre et al., 2018; it is increasingly clear that individuals with ablated gills have lower metabolic rate and maximal metabolic rate may be depressed in larger individuals in warmer (Pauly, 2010;Rubalcaba et al., 2020). Assuming that (gill) oxygen limitation is stronger in larger species than smaller ones, oxygen uptake might not be limited in smaller species where other factors might be shaping r max . Thus, the positive effect of temperature on r max in smaller species might be a direct result of the increased speed of chemical and metabolic reactions as hypothesized by metabolic theory Savage et al., 2004). In larger species, the positive effect of temperature on metabolic reactions could be completely overridden by the oxygen limitation that gills pose as organisms increase in volume (Pauly & Cheung, 2018), resulting in a negative relationship between temperature and r max . Nonetheless, some species might be behaviorally evading these temperature (and oxygen) constraints. For example, the Whale Shark (Rhincodon typus) is the largest extant fish and lives in warm waters throughout the world's oceans, yet it performs deep dives (>1000 m) into cold waters, which are hypothesized to aid in thermoregulation to dissipate excess heat (Thums et al., 2012), and our finding leads us to speculate that they perhaps also access waters with higher oxygen concentrations. While the effect of oxygen concentrations across depth on r max was not explored in this manuscript, it might further explain the low r max estimates seen among deep-sea species in line with the expectations from the GOLT and other general theories of oxygen limitation Pauly, 2021;Pörtner et al., 2017;Rubalcaba et al., 2020).
The existence of a depth dimension to r max , albeit complex, is worth exploring further. Many potential mechanisms could drive the decrease in r max among chondrichthyans with increasing depth, all of which relate to the amount of energy available for metabolic processes: temperature, light, and consequently primary production decrease below the photic zone (Gage & Tyler, 1991;Jahnke, 1996); metabolic capacity and requirements decrease with depth as dissolved oxygen and activity levels decrease (Childress et al., 1990;; animal biomass, and hence food availability, also decreases with increasing depth (Rex et al., 2006); even the unique physiology of sharks and rays can increase the energetic cost of living in the deep, and consequently reduce the energy available for production (Treberg & Speers-Roesch, 2016 TA B L E 2 Comparison of log(r max ) models using standard and corrected Akaike information criteria (AICc), number of parameters (n), negative log-likelihood (−LL), adjusted R 2 , and Akaike weights.

F I G U R E 2
Coefficient plots for the four models of log(r max ) with lowest AICc values. Lighter colors indicate models with decreasing support based on ΔAICc. Error bars show the 95% confidence intervals, and effect sizes were considered significant when confidence intervals do not overlap zero. Shaded areas show the expected effect sizes for body mass (−0.33 to −0.25) and temperature (−1.0 to −0.6) based on metabolic scaling theory. , rather than lower food availability or other factors.
This hypothesis posits that deep-sea species have lower activity levels as they interact with predators and prey on much smaller spatial scales, or less frequently, than shallow water ones, as the distance at which predators and prey interact is reduced as light levels drop, resulting in lower basal metabolic rates and hence lower productivities , and also shallower mass scaling of metabolic rate in deeper waters (Rubalcaba et al., 2020). A decrease in r max with increasing depth could also arise from the unique physiological

F I G U R E 3
Relationship between maximum weight and maximum intrinsic rate of population increase r max , in log space, for 63 chondrichthyan species. Median depth and temperature for each species are shown by the point size and color, respectively. Median temperatures are corrected for species, which have body temperatures that are higher than their surroundings. Fitted lines show predicted relationships based on the top-ranked model. (a) Predicted allometric changes of r max across median depths (10, 500, 1000 m) but constant median temperature (6°C), and (b) predicted allometric changes of r max for three different median temperatures (6, 10, 20°C) but constant median depth (10 m). energetic gradients discussed above can potentially affect the productivity and population growth rates of species found at a different depth, even after the temperature is accounted for, and as such, depth can be thought of as an added axis of variation in the slow-fast lifehistory continuum (Sibly & Brown, 2007).
The strong phylogenetic signal in the residuals of r max suggests that other aspects of biology shape their maximum population growth rate, which is likely evolutionarily conserved traits. This likely stems, in part, from estimates of r max being strongly influenced by reproductive output (Pardo et al., 2018), which is strongly conserved among closely related species (Dulvy & Reynolds, 1997). This strong phylogenetic signal opens up the door for predictive modeling of r max based on phylogenetic relationships.
While our study is a first attempt at understanding the relationship between r max , size, and the environment, there is a considerable amount of unexplained variation in our models. Considering the caveats when interpreting the findings of this study might enable future studies to improve model fits. First, there is a considerable amount of collinearity between depth and temperature (Pearson's r = .64). This is somewhat unavoidable as warmer waters are almost chiefly found in shallow depths, while deep waters are consistently cold. While the diagnostic tests we performed suggest that our results are robust to this collinearity, it is important to bear this in mind, particularly if our framework is to be adapted in future for predictive purposes. Second, our study is limited by the choice of hypotheses being examined. While we completed nine different hypotheses using an information-theoretic approach, it is possible that other hypotheses not considered in this study are better at explaining the patterns in the data. Including other environmental variables mentioned above, such as dissolved oxygen or net primary production, could result in models that are better supported in the data than the models compared in this study. Concomitantly, we only explored simple interactions from a linear modeling framework; other nonlinear frameworks for modeling might be better supported by the data. Third, we used maximum weight as our metric of mass as it is the most readily available trait among data-poor sharks and rays; however, it is also among the most variable and our choice of this life-history trait might contribute to the low R 2 of our models. Weight at maturity is a likely better parameter to use, albeit more difficult to estimate, as it is more precise than estimates of maximum weight, and also aligns better with expectations of metabolic scaling  and life-history theories (Charnov et al., 2007). Lastly, our method of accounting for mesothermy (i.e., incrementing environmental temperature by a fixed amount for all mesotherms) might not be the most adequate; there are likely large interspecific differences in the degree of mesothermy, and also how the temperature differential varies across a range of environmental temperatures.
Our study explored the relationships between body mass, temperature, and depth with the maximum intrinsic rate of increase in sharks and rays and paves the way for further investigation of the mechanisms behind these relationships. For example, investigating how the mass scaling relationship in chondrichthyans from solely shallow environments TA B L E 4 Differences in corrected Akaike Information Criteria (ΔAICc) for the models run with 20 different iterations of the phylogenetic tree published by Stein et al. (2018). log  varies with temperature might confirm that the apparent "flattening" of the mass scaling relationship of productivity at depth was indeed due to temperature. Taking into account additional spatial datasets that provide information on nutrient availability in the deep sea, such as sedimentation rates or organic carbon burial rates (Jahnke, 1996), could provide a better understanding of the degree to which food availability limits productivity among deep-sea sharks and rays. Similarly, incorporating data on dissolved oxygen content would elucidate whether the observed temperature effects, particularly among medium to largesized species, are a result of oxygen limitation. Alternatively, exploring in more detail whether the effect of depth on r max levels off after a certain threshold (as it does for basal metabolic rate in pelagic fishes, supporting the visual-interaction hypothesis) or whether it is continuous would help elucidate the mechanism driving the observed relationship with depth. Identifying the mechanisms behind these observed relationships among body size, temperature, and depth will eventually provide a much clearer understanding of how the productivity, and therefore the vulnerability of species, varies across the world's oceans.

ACK N OWLED G M ENTS
We are grateful to Jennifer Bigman, Philina English, Sarah Gravel,

CO N FLI C T O F I NTE R E S T
None.